Glacier extinction homogenizes functional diversity
Supplementary material - R script
Khelidj N, Cerabolini BEL, Caccianiga M, Losapio G
Introduction
This Supplemetary Material reports the code to reproduce the analyses and figures carried out in the manuscript entitled “Glacier extinction homogenizes functional diversity” by Khelidj, Cerabolini, Caccianiga & Losapio, 2023, Journal of Vegetation Science.
The scope of this tutorial is to increase its reproducibility, clarity, transparency and dissemination. Data analysis was done using R version 4.1.2. This document was compiled with the `rmarkdown package, version Version 1.2.5033. This tutorial is licensed under CC BY-NC-ND 4.0, which means you are free to share, copy and redistribute this tutorial in any medium or format under the terms of attribution of appropriate credit, non-commercial purposes and no derivatives.
The citation is:
Khelidj R, Cerabolini BEL, Caccianiga M, Losapio G. 2023. Glacier extinction homogenizes functional diversity. Journal of Vegetation Science.
Software preparation
Install R (https://www.r-project.org) (R Core Team 2022) if you do not have it already. Then, install and load the following packages.
library(readxl)
library(FD)
library(ggplot2)
library(ggpubr)
library(lme4)
library(lmerTest)
library(car)
library(effects)
library(effectsize)
library(parameters)
library(merDeriv)
library(mgcv)
library(adiv)
library(maditr)Data preparation
Download and import the following data tables
## trait data
traitsdb = read.csv("TraitsDB.csv")
rownames(traitsdb) = traitsdb$species
traitsdb$X = NULL
# species occurrence data
# Amola glacier
amola_matb = read.csv('amola_matb.csv')
amola_env = read.csv('amola_env.csv')
# Cedec glacier
cedec_matb = read.csv("cedec_matb.csv")
cedec_env = read.csv("cedec_envb.csv")
# Rutor glacier
rutor_matb = read.csv("rutor_matb.csv")
rutor_matb = read.csv("rutor_env.csv")
# Trobio glacier
trobio_matb = read.csv("trobio_matb.csv")
trobio_matb = read.csv("trobio_env.csv")
###
matb = list(amola_matb,cedec_matb,rutor_matb,trobio_matb)Transform glacier retreat data
amola_env$lyears = log(amola_env$years+1)
cedec_env$lyears = log(cedec_env$years+1)
rutor_env$lyears = log(rutor_env$years+1)
trobio_env$lyears = log(trobio_env$years+1)
###
matenv = list(amola_env,cedec_env,rutor_env,trobio_env)Calculate CV
### coefficient of variation
traitsdb$LNC.se = ifelse(traitsdb$LNC.se==0, NA, traitsdb$LNC.se)
traitsdb$LCC.se = ifelse(traitsdb$LCC.se==0, NA, traitsdb$LCC.se)
traitsdb$CN.se = ifelse(traitsdb$CN.se==0, NA, traitsdb$CN.se)
traitsdb$LA.se = ifelse(traitsdb$LA.se==0, NA, traitsdb$LA.se)
traitsdb$LDW.se = ifelse(traitsdb$LDW.se==0, NA, traitsdb$LDW.se)
traitsdb$LFW.se = ifelse(traitsdb$LFW.se==0, NA, traitsdb$LFW.se)
traitsdb$SLA.se = ifelse(traitsdb$SLA.se==0, NA, traitsdb$SLA.se)
traitsdb$LDMC.se = ifelse(traitsdb$LDMC.se==0, NA, traitsdb$LDMC.se)
traitsdb$CAN.se = ifelse(traitsdb$CAN.se==0, NA, traitsdb$CAN.se)
#
traitsdb$lnc_cv = 1/traitsdb$LNC.med * sqrt(traitsdb$LNC.n) * traitsdb$LNC.se
traitsdb$lcc_cv = 1/traitsdb$LCC.med * sqrt(traitsdb$LCC.n) * traitsdb$LCC.se
traitsdb$cn_cv = 1/traitsdb$CN.med * sqrt(traitsdb$CN.n) * traitsdb$CN.se
traitsdb$la_cv = 1/traitsdb$LA.med * sqrt(traitsdb$LA.n) * traitsdb$LA.se
traitsdb$ldw_cv = 1/traitsdb$LDW.med * sqrt(traitsdb$LDW.n) * traitsdb$LDW.se
traitsdb$ldw_cv[which(traitsdb$ldw_cv==Inf)] = NA
traitsdb$lfw_cv = 1/traitsdb$LFW.med * sqrt(traitsdb$LFW.n) * traitsdb$LFW.se
traitsdb$lfw_cv[which(traitsdb$lfw_cv==Inf)] = NA
traitsdb$sla_cv = 1/traitsdb$SLA.med * sqrt(traitsdb$SLA.n) * traitsdb$SLA.se
traitsdb$ldmc.cv = 1/traitsdb$LDMC.med * sqrt(traitsdb$LDMC.n) * traitsdb$LDMC.se
traitsdb$can_cv = 1/traitsdb$CAN.med * sqrt(traitsdb$CAN.n) * traitsdb$CAN.seSynthesis indices
We calculate and syntheses functional diversity across plant traits, considering both trait average and trait variation.
Functional evenness and functional dispersion
We do so my means of distance-based functional diversity indices that
address functional evenness (FEve) and functional dispersion (FDis). We
use the R function dbFD in FD package.
for(i in 1:4){
print(i)
### mean ###
traits_mn = traitsdb[which(rownames(traitsdb)%in%colnames(matb[[i]])),
c(2,3,4,5,11,14,17,20,23,26,29,32,38)]
##### observed trends in average traits
fcomp = dbFD(traits_mn, matb[[i]], corr = 'none')
if(i==1){
alltr_mni = data.frame(cbind(matenv[[i]]$lyears, fcomp$FEve, fcomp$FDis,
rep(i,nrow(matenv[[i]]))))
}
if(i>1){
m0 = data.frame(cbind(matenv[[i]]$lyears, fcomp$FEve, fcomp$FDis,
rep(i,nrow(matenv[[i]]))))
alltr_mni = data.frame(rbind(alltr_mni,m0))
}
### cv ###
traits_cv = traitsdb[which(rownames(traitsdb)%in%colnames(matb[[i]])),
41:49]
##### observed trends in trait cv
fcomp = dbFD(traits_cv, matb[[i]], corr = 'none')
if(i==1){
alltr_cvi = data.frame(cbind(matenv[[i]]$lyears, fcomp$FEve, fcomp$FDis,
rep(i,nrow(matenv[[i]]))))
}
if(i>1){
m0 = data.frame(cbind(matenv[[i]]$lyears, fcomp$FEve, fcomp$FDis,
rep(i,nrow(matenv[[i]]))))
alltr_cvi = data.frame(rbind(alltr_cvi,m0))
}
}
#
colnames(alltr_mni)[1] = 'years'
colnames(alltr_mni)[2] = 'feve'
colnames(alltr_mni)[3] = 'fdis'
colnames(alltr_mni)[4] = 'site'
alltr_mni$site = as.factor(alltr_mni$site)
alltr_mni$sprich = 0
colnames(alltr_cvi)[1] = 'years'
colnames(alltr_cvi)[2] = 'feve'
colnames(alltr_cvi)[3] = 'fdis'
colnames(alltr_cvi)[4] = 'site'
alltr_cvi$site = as.factor(alltr_cvi$site)Evaluate relationship with species richness
alltr_cvi$sprich = 0
### relationship with alpha diversity ####
for(i in 1:4){
print(i)
alltr_mni$sprich[which(alltr_mni$site==i)] = rowSums(matb[[i]])
alltr_cvi$sprich[which(alltr_cvi$site==i)] = rowSums(matb[[i]])
}Run statistical analyses
We then look how functional evenness and functional dispersion change with glacier retreat.
Functional evenness of trait average
#### stat model ###
# mean feve
mod_feve_mni = lmer (feve ~ poly(years,2)*sprich + (1|site),
data = alltr_mni)
anova(mod_feve_mni)
cohens_f(mod_feve_mni, partial = TRUE)
modpars = model_parameters(mod_feve_mni, effects = "fixed", ci_method = "satterthwaite")
effsz = t_to_r(c(modpars$t[1],modpars$t[2],modpars$t[3],
modpars$t[4],modpars$t[5],modpars$t[6]),
df_error = c(modpars$df[1],modpars$df[2],modpars$df[3],
modpars$df[4],modpars$df[5],modpars$df[6]))Print model output and effect size
options(width = 100)
print(modpars)## # Fixed Effects
##
## Parameter | Coefficient | SE | 95% CI | t | df | p
## -----------------------------------------------------------------------------------------------
## (Intercept) | 0.29 | 0.03 | [ 0.21, 0.37] | 8.49 | 7.74 | < .001
## years [1st degree] | -1.47 | 0.24 | [-1.95, -0.99] | -6.10 | 139.96 | < .001
## years [2nd degree] | 1.47 | 0.18 | [ 1.12, 1.82] | 8.35 | 141.39 | < .001
## sprich | -9.79e-03 | 1.51e-03 | [-0.01, -0.01] | -6.48 | 125.24 | < .001
## years [1st degree] * sprich | 0.10 | 0.02 | [ 0.06, 0.13] | 5.68 | 140.14 | < .001
## years [2nd degree] * sprich | -0.09 | 0.02 | [-0.12, -0.06] | -5.93 | 140.89 | < .001##
## Uncertainty intervals (equal-tailed) and p-values (two-tailed) computed using a Wald
## t-distribution with Satterthwaite approximation.print(effsz)## r | 95% CI
## ----------------------
## 0.95 | [ 0.81, 0.98]
## -0.46 | [-0.57, -0.32]
## 0.57 | [ 0.46, 0.66]
## -0.50 | [-0.61, -0.36]
## 0.43 | [ 0.29, 0.55]
## -0.45 | [-0.56, -0.31]Functional evenness of trait variation
# cv feve
mod_feve_cvi = lmer (feve ~ poly(years,2)*sprich + (1|site),
data = alltr_cvi)
anova(mod_feve_cvi)
cohens_f(mod_feve_cvi, partial = TRUE)
modpars = model_parameters(mod_feve_cvi,
effects = "fixed",
ci_method = "satterthwaite")
effsz = t_to_r(c(modpars$t[1],modpars$t[2],modpars$t[3],
modpars$t[4],modpars$t[5],modpars$t[6]),
df_error = c(modpars$df[1],modpars$df[2],modpars$df[3],
modpars$df[4],modpars$df[5],modpars$df[6]))Print model output and effect size
print(modpars)## # Fixed Effects
##
## Parameter | Coefficient | SE | 95% CI | t | df | p
## -----------------------------------------------------------------------------------------------
## (Intercept) | 0.29 | 0.03 | [ 0.21, 0.37] | 8.49 | 7.74 | < .001
## years [1st degree] | -1.47 | 0.24 | [-1.95, -0.99] | -6.10 | 139.96 | < .001
## years [2nd degree] | 1.47 | 0.18 | [ 1.12, 1.82] | 8.35 | 141.39 | < .001
## sprich | -9.79e-03 | 1.51e-03 | [-0.01, -0.01] | -6.48 | 125.24 | < .001
## years [1st degree] * sprich | 0.10 | 0.02 | [ 0.06, 0.13] | 5.68 | 140.14 | < .001
## years [2nd degree] * sprich | -0.09 | 0.02 | [-0.12, -0.06] | -5.93 | 140.89 | < .001##
## Uncertainty intervals (equal-tailed) and p-values (two-tailed) computed using a Wald
## t-distribution with Satterthwaite approximation.print(effsz)## r | 95% CI
## ----------------------
## 0.95 | [ 0.81, 0.98]
## -0.46 | [-0.57, -0.32]
## 0.57 | [ 0.46, 0.66]
## -0.50 | [-0.61, -0.36]
## 0.43 | [ 0.29, 0.55]
## -0.45 | [-0.56, -0.31]Functional dispersion of trait average
# mean fdis
mod_fdis_mni = lmer (fdis ~ poly(years,2)*sprich + (1|site),
data = alltr_mni)
anova(mod_fdis_mni)
cohens_f(mod_fdis_mni, partial = TRUE)
modpars = model_parameters(mod_fdis_mni,
effects = "fixed",
ci_method = "satterthwaite")
effsz = t_to_r(c(modpars$t[1],modpars$t[2],modpars$t[3],
modpars$t[4],modpars$t[5],modpars$t[6]),
df_error = c(modpars$df[1],modpars$df[2],modpars$df[3],
modpars$df[4],modpars$df[5],modpars$df[6]))Print model output and effect size
options(width = 100)
print(modpars)## # Fixed Effects
##
## Parameter | Coefficient | SE | 95% CI | t | df | p
## -----------------------------------------------------------------------------------------------
## (Intercept) | 0.29 | 0.03 | [ 0.21, 0.37] | 8.49 | 7.74 | < .001
## years [1st degree] | -1.47 | 0.24 | [-1.95, -0.99] | -6.10 | 139.96 | < .001
## years [2nd degree] | 1.47 | 0.18 | [ 1.12, 1.82] | 8.35 | 141.39 | < .001
## sprich | -9.79e-03 | 1.51e-03 | [-0.01, -0.01] | -6.48 | 125.24 | < .001
## years [1st degree] * sprich | 0.10 | 0.02 | [ 0.06, 0.13] | 5.68 | 140.14 | < .001
## years [2nd degree] * sprich | -0.09 | 0.02 | [-0.12, -0.06] | -5.93 | 140.89 | < .001##
## Uncertainty intervals (equal-tailed) and p-values (two-tailed) computed using a Wald
## t-distribution with Satterthwaite approximation.print(effsz)## r | 95% CI
## ----------------------
## 0.95 | [ 0.81, 0.98]
## -0.46 | [-0.57, -0.32]
## 0.57 | [ 0.46, 0.66]
## -0.50 | [-0.61, -0.36]
## 0.43 | [ 0.29, 0.55]
## -0.45 | [-0.56, -0.31]Functional dispersion of trait variation
# cv fdis
mod_fdis_cvi = lmer (fdis ~ poly(years,2)*sprich + (1|site),
data = alltr_cvi)
anova(mod_fdis_cvi)
cohens_f(mod_fdis_cvi, partial = TRUE)
modpars = model_parameters(mod_fdis_cvi,
effects = "fixed",
ci_method = "satterthwaite")
effsz = t_to_r(c(modpars$t[1],modpars$t[2],modpars$t[3],
modpars$t[4],modpars$t[5],modpars$t[6]),
df_error = c(modpars$df[1],modpars$df[2],modpars$df[3],
modpars$df[4],modpars$df[5],modpars$df[6]))Print model output and effect size
print(modpars)## # Fixed Effects
##
## Parameter | Coefficient | SE | 95% CI | t | df | p
## -----------------------------------------------------------------------------------------------
## (Intercept) | 0.29 | 0.03 | [ 0.21, 0.37] | 8.49 | 7.74 | < .001
## years [1st degree] | -1.47 | 0.24 | [-1.95, -0.99] | -6.10 | 139.96 | < .001
## years [2nd degree] | 1.47 | 0.18 | [ 1.12, 1.82] | 8.35 | 141.39 | < .001
## sprich | -9.79e-03 | 1.51e-03 | [-0.01, -0.01] | -6.48 | 125.24 | < .001
## years [1st degree] * sprich | 0.10 | 0.02 | [ 0.06, 0.13] | 5.68 | 140.14 | < .001
## years [2nd degree] * sprich | -0.09 | 0.02 | [-0.12, -0.06] | -5.93 | 140.89 | < .001##
## Uncertainty intervals (equal-tailed) and p-values (two-tailed) computed using a Wald
## t-distribution with Satterthwaite approximation.print(effsz)## r | 95% CI
## ----------------------
## 0.95 | [ 0.81, 0.98]
## -0.46 | [-0.57, -0.32]
## 0.57 | [ 0.46, 0.66]
## -0.50 | [-0.61, -0.36]
## 0.43 | [ 0.29, 0.55]
## -0.45 | [-0.56, -0.31]Make figures
##### plots #####
plot(Effect(c("years"), mod_feve_mni))
plot(Effect(c("years"), mod_fdis_mni))
plot(Effect(c("years"), mod_feve_cvi))
plot(Effect(c("years"), mod_fdis_cvi))
plot(Effect(c("sprich"), mod_feve_mni))
plot(Effect(c("sprich"), mod_fdis_mni))
plot(Effect(c("sprich"), mod_feve_cvi))
plot(Effect(c("sprich"), mod_fdis_cvi))
plot(Effect(c("sprich", "years"), mod_feve_mni))
plot(Effect(c("sprich", "years"), mod_fdis_mni))
plot(Effect(c("sprich", "years"), mod_feve_cvi))
plot(Effect(c("sprich", "years"), mod_fdis_cvi))
### gg plots
a1 = ggplot(alltr_mni, aes(x = years,
y = feve,
colour = site)) +
geom_jitter(width = 0.2) +
geom_smooth(method = lm, formula = y ~ poly(x,2)) +
theme_classic()
a2 = ggplot(alltr_mni, aes(x = years,
y = fdis,
colour = site)) +
geom_jitter(width = 0.2) +
geom_smooth(method = lm, formula = y ~ poly(x,2)) +
theme_classic()
a3 = ggplot(alltr_cvi, aes(x = years,
y = feve,
colour = site)) +
geom_jitter(width = 0.2) +
geom_smooth(method = lm, formula = y ~ poly(x,2)) +
theme_classic()
a4 = ggplot(alltr_cvi, aes(x = years,
y = fdis,
colour = site)) +
geom_jitter(width = 0.2) +
geom_smooth(method = lm, formula = y ~ poly(x,2)) +
theme_classic()
ggarrange(a1, a3, a2, a4,
ncol = 2, nrow = 2)## Warning: Removed 10 rows containing non-finite values (stat_smooth).## Warning: Removed 10 rows containing missing values (geom_point).## Warning: Removed 10 rows containing non-finite values (stat_smooth).## Warning: Removed 10 rows containing missing values (geom_point).
### sprich - fdiv
##### plots #####
a1 = ggplot(alltr_mni, aes(x = sprich,
y = feve,
colour = site)) +
geom_jitter(width = 0.2) +
geom_smooth(method = lm, formula = y ~ poly(x,2)) +
theme_classic()
a2 = ggplot(alltr_mni, aes(x = sprich,
y = fdis,
colour = site)) +
geom_jitter(width = 0.2) +
geom_smooth(method = lm, formula = y ~ poly(x,2)) +
theme_classic()
a3 = ggplot(alltr_cvi, aes(x = sprich,
y = feve,
colour = site)) +
geom_jitter(width = 0.2) +
geom_smooth(method = lm, formula = y ~ poly(x,2)) +
theme_classic()
a4 = ggplot(alltr_cvi, aes(x = sprich,
y = fdis,
colour = site)) +
geom_jitter(width = 0.2) +
geom_smooth(method = lm, formula = y ~ poly(x,2)) +
theme_classic()
ggarrange(a1, a3, a2, a4,
ncol = 2, nrow = 2)## Warning: Removed 10 rows containing non-finite values (stat_smooth).
## Removed 10 rows containing missing values (geom_point).## Warning: Removed 10 rows containing non-finite values (stat_smooth).## Warning: Removed 10 rows containing missing values (geom_point).
Functional redundancy and uniquenness
We then consider functional redundancy and uniquenness for each glacier foreland site.
#### functional redundancy ########
## amola
amola_distm = gowdis(traits_mn[which(rownames(traits_mn)%in%colnames(amola_matb)),])
amola_uniq = uniqueness(amola_matb[,which(colnames(amola_matb)%in%rownames(traits_mn))],
amola_distm,
abundance = F)
str(amola_uniq)
amola_env = cbind(amola_env,amola_uniq$red[,c(4,6)])
## cedec
cedec_distm = gowdis(traits_mn[which(rownames(traits_mn)%in%colnames(cedec_matb)),])
cedec_uniq = uniqueness(cedec_matb[,which(colnames(cedec_matb)%in%rownames(traits_mn))],
cedec_distm,
abundance = F)
cedec_uniq$Kbar
cedec_env = cbind(cedec_env,cedec_uniq$red[,c(4,6)])
## rutor
rutor_distm = gowdis(traits_mn[which(rownames(traits_mn)%in%colnames(rutor_matb)),])
rutor_uniq = uniqueness(rutor_matb[,which(colnames(rutor_matb)%in%rownames(traits_mn))],
rutor_distm,
abundance = F)
rutor_uniq$Kbar
rutor_env = cbind(rutor_env,rutor_uniq$red[,c(4,6)])
## trobio
trobio_distm = gowdis(traits_mn[which(rownames(traits_mn)%in%colnames(trobio_matb)),])
trobio_uniq = uniqueness(trobio_matb[,which(colnames(trobio_matb)%in%rownames(traits_mn))],
trobio_distm,
abundance = F)
trobio_uniq$Kbar
trobio_env = cbind(trobio_env,trobio_uniq$red[,c(4,6)])
### pool data
alltr_mni$ustar = c(amola_env$Ustar, cedec_env$Ustar, rutor_env$Ustar, trobio_env$Ustar)Statistical analysis
mod_ustar = lmer (ustar ~ poly(years,2)*sprich + (1|site),
data = alltr_mni)
anova(mod_ustar)
cohens_f(mod_ustar, partial = TRUE)
modpars = model_parameters(mod_ustar, effects = "fixed", ci_method = "satterthwaite")
effsz = t_to_r(c(modpars$t[1],modpars$t[2],modpars$t[3],
modpars$t[4],modpars$t[5],modpars$t[6]),
df_error = c(modpars$df[1],modpars$df[2],modpars$df[3],
modpars$df[4],modpars$df[5],modpars$df[6]))Print model output
print(modpars)## # Fixed Effects
##
## Parameter | Coefficient | SE | 95% CI | t | df | p
## -----------------------------------------------------------------------------------------------
## (Intercept) | 0.29 | 0.03 | [ 0.21, 0.37] | 8.49 | 7.74 | < .001
## years [1st degree] | -1.47 | 0.24 | [-1.95, -0.99] | -6.10 | 139.96 | < .001
## years [2nd degree] | 1.47 | 0.18 | [ 1.12, 1.82] | 8.35 | 141.39 | < .001
## sprich | -9.79e-03 | 1.51e-03 | [-0.01, -0.01] | -6.48 | 125.24 | < .001
## years [1st degree] * sprich | 0.10 | 0.02 | [ 0.06, 0.13] | 5.68 | 140.14 | < .001
## years [2nd degree] * sprich | -0.09 | 0.02 | [-0.12, -0.06] | -5.93 | 140.89 | < .001##
## Uncertainty intervals (equal-tailed) and p-values (two-tailed) computed using a Wald
## t-distribution with Satterthwaite approximation.print(effsz)## r | 95% CI
## ----------------------
## 0.95 | [ 0.81, 0.98]
## -0.46 | [-0.57, -0.32]
## 0.57 | [ 0.46, 0.66]
## -0.50 | [-0.61, -0.36]
## 0.43 | [ 0.29, 0.55]
## -0.45 | [-0.56, -0.31]Figures
### Make figures
ggplot(amola_env, aes(x = lyears,
y = Ustar)) +
geom_jitter(width = 0.2) +
geom_smooth(method = lm, formula = y ~ poly(x,2)) +
theme_classic()## Warning: Removed 1 rows containing non-finite values (stat_smooth).## Warning: Removed 1 rows containing missing values (geom_point).
ggplot(cedec_env, aes(x = lyears,
y = Ustar)) +
geom_jitter(width = 0.2) +
geom_smooth(method = lm, formula = y ~ poly(x,2)) +
theme_classic()
ggplot(rutor_env, aes(x = lyears,
y = Ustar)) +
geom_jitter(width = 0.2) +
geom_smooth(method = lm, formula = y ~ poly(x,2)) +
theme_classic()
ggplot(trobio_env, aes(x = lyears,
y = Ustar)) +
geom_jitter(width = 0.2) +
geom_smooth(method = lm, formula = y ~ poly(x,2)) +
theme_classic()
ggplot(alltr_mni, aes(x = years,
y = ustar,
colour = site)) +
geom_jitter(width = 0.2) +
geom_smooth(method = lm, formula = y ~ poly(x,2)) +
theme_classic()## Warning: Removed 1 rows containing non-finite values (stat_smooth).
## Removed 1 rows containing missing values (geom_point).
Functional diversity at the trait level
We then look at trait average and trait variation for each single trait separately.
for(i in 1:4){
print(i)
### mean ###
traits_mn = traitsdb[which(rownames(traitsdb)%in%colnames(matb[[i]])),
c(3,4,5,11,14,17,20,23,26,29,32,38)]
##### observed trends in average traits
fcomp = functcomp(traits_mn, matb[[i]])
if(i==1){
alltr_mn = cbind(matenv[[i]]$lyears, fcomp, rep(i,nrow(matenv[[i]])))
}
if(i>1){
m0 = cbind(matenv[[i]]$lyears, fcomp, rep(i,nrow(matenv[[i]])))
alltr_mn = rbind(alltr_mn,m0)
}
### cv ###
traits_cv = traitsdb[which(rownames(traitsdb)%in%colnames(matb[[i]])),
41:49]
##### observed trends in trait cv
fcomp = functcomp(traits_cv, matb[[i]])
if(i==1){
alltr_cv = cbind(matenv[[i]]$lyears, fcomp, rep(i,nrow(matenv[[i]])))
}
if(i>1){
m0 = cbind(matenv[[i]]$lyears, fcomp, rep(i,nrow(matenv[[i]])))
alltr_cv = rbind(alltr_cv,m0)
}
}
#
colnames(alltr_mn)[1] = 'years'
colnames(alltr_mn)[2] = 'latspr'
colnames(alltr_mn)[3] = 'flost'
colnames(alltr_mn)[4] = 'floper'
colnames(alltr_mn)[5] = 'lnc'
colnames(alltr_mn)[6] = 'lcc'
colnames(alltr_mn)[7] = 'cn'
colnames(alltr_mn)[8] = 'la'
colnames(alltr_mn)[9] = 'ldw'
colnames(alltr_mn)[10] = 'lfw'
colnames(alltr_mn)[11] = 'sla'
colnames(alltr_mn)[12] = 'ldmc'
colnames(alltr_mn)[13] = 'can'
colnames(alltr_mn)[14] = 'site'
alltr_mn$site = as.factor(alltr_mn$site)
#
colnames(alltr_cv)[1] = 'years'
colnames(alltr_cv)[2] = 'lnc'
colnames(alltr_cv)[3] = 'lcc'
colnames(alltr_cv)[4] = 'cn'
colnames(alltr_cv)[5] = 'la'
colnames(alltr_cv)[6] = 'ldw'
colnames(alltr_cv)[7] = 'lfw'
colnames(alltr_cv)[8] = 'sla'
colnames(alltr_cv)[9] = 'ldmc'
colnames(alltr_cv)[10] = 'can'
colnames(alltr_cv)[11] = 'site'
alltr_cv$site = as.factor(alltr_cv$site)Statistical analysis
## mixed models & effect size####
## mean
mod.alltr.mn = rep(list(NA),12)
st.alltr.mn = data.frame(trait = gl(12,2,24, labels=colnames(alltr_mn)[2:13]),
estimate = NA,
stderr = NA,
bcil = NA,
bciu = NA,
df = NA,
tval = NA,
pval = NA,
choenr = NA,
rcil = NA,
rciu = NA)
for(i in 1:12){
print(i)
flamod = as.formula(paste(paste(colnames(alltr_mn)[i+1], "~"),
paste('poly(years,2)', '(1|site)', sep="+")))
mod.alltr.mn[[i]] = lmer(flamod, data = alltr_mn)
modpars = model_parameters(mod.alltr.mn[[i]],
effects = "fixed",
ci_method = "satterthwaite")
st.alltr.mn[((i-1)*2+1),2:8] = as.numeric(modpars[2,c(2:9)])[-2]
st.alltr.mn[((i-1)*2+2),2:8] = as.numeric(modpars[3,c(2:9)])[-2]
effsz = t_to_r(c(modpars$t[2],modpars$t[3]),
df_error = c(modpars$df[2],modpars$df[3]))
st.alltr.mn[((i-1)*2+1),9:11] = as.numeric(effsz[1,])[-2]
st.alltr.mn[((i-1)*2+2),9:11] = as.numeric(effsz[2,])[-2]
}
## cv
mod.alltr.cv = rep(list(NA),9)
st.alltr.cv = data.frame(trait = gl(9,2,18, labels=colnames(alltr_cv)[2:10]),
estimate = NA,
stderr = NA,
bcil = NA,
bciu = NA,
df = NA,
tval = NA,
pval = NA,
choenr = NA,
rcil = NA,
rciu = NA)
for(i in 1:9){
print(i)
if(i!=1 & i!=6){
flamod = as.formula(paste(paste(colnames(alltr_cv)[i+1], "~"),
paste('poly(years,2)', '(1|site)', sep="+")))
mod.alltr.cv[[i]] = lmer(flamod, data = alltr_cv)
modpars = model_parameters(mod.alltr.cv[[i]],
effects = "fixed",
ci_method = "satterthwaite")
st.alltr.cv[((i-1)*2+1),2:8] = as.numeric(modpars[2,c(2:9)])[-2]
st.alltr.cv[((i-1)*2+2),2:8] = as.numeric(modpars[3,c(2:9)])[-2]
effsz = t_to_r(c(modpars$t[2],modpars$t[3]),
df_error = c(modpars$df[2],modpars$df[3]))
st.alltr.cv[((i-1)*2+1),9:11] = as.numeric(effsz[1,])[-2]
st.alltr.cv[((i-1)*2+2),9:11] = as.numeric(effsz[2,])[-2]}
if(i==1 | i==6){
flamod = as.formula(paste(paste(colnames(alltr_cv)[i+1], "~"),
paste('years', '(1|site)', sep="+")))
mod.alltr.cv[[i]] = lmer(flamod, data = alltr_cv)
modpars = model_parameters(mod.alltr.cv[[i]], effects = "fixed", ci_method = "satterthwaite")
st.alltr.cv[((i-1)*2+1),2:8] = as.numeric(modpars[2,c(2:9)])[-2]
effsz = t_to_r(modpars$t[2],modpars$df[2])
st.alltr.cv[((i-1)*2+1),9:11] = as.numeric(effsz[1,])[-2]}
}Print model outcome
## summary table
st.alltr.mn[,-1] = round(st.alltr.mn[,-1],3)
st.alltr.cv[,-1] = round(st.alltr.cv[,-1],3)
print(st.alltr.mn)## trait estimate stderr bcil bciu df tval pval choenr rcil rciu
## 1 latspr -1.004 0.95 -1.540 -0.467 -3.696 144.220 0.000 -0.294 -0.429 -0.139
## 2 latspr 0.560 0.95 -0.018 1.137 1.915 145.913 0.057 0.157 -0.005 0.306
## 3 flost -0.825 0.95 -1.149 -0.500 -5.021 143.573 0.000 -0.386 -0.507 -0.241
## 4 flost 0.005 0.95 -0.347 0.356 0.026 145.153 0.979 0.002 -0.158 0.163
## 5 floper 0.537 0.95 0.289 0.784 4.279 143.431 0.000 0.336 0.185 0.465
## 6 floper -0.158 0.95 -0.426 0.111 -1.159 144.784 0.248 -0.096 -0.251 0.067
## 7 lnc 0.156 0.95 -0.387 0.700 0.569 144.166 0.571 0.047 -0.115 0.206
## 8 lnc -1.284 0.95 -1.869 -0.699 -4.338 145.952 0.000 -0.338 -0.465 -0.188
## 9 lcc 9.129 0.95 7.709 10.550 12.704 144.381 0.000 0.726 0.649 0.783
## 10 lcc 3.305 0.95 1.778 4.831 4.278 145.664 0.000 0.334 0.183 0.462
## 11 cn 7.331 0.95 1.987 12.675 2.712 143.848 0.008 0.221 0.060 0.364
## 12 cn 14.132 0.95 8.361 19.903 4.840 145.824 0.000 0.372 0.226 0.494
## 13 la 561.067 0.95 237.151 884.983 3.424 143.920 0.001 0.274 0.117 0.412
## 14 la 54.999 0.95 -294.804 404.802 0.311 145.821 0.756 0.026 -0.135 0.185
## 15 ldw 37.559 0.95 19.369 55.750 4.081 143.901 0.000 0.322 0.169 0.453
## 16 ldw 2.393 0.95 -17.255 22.042 0.241 145.786 0.810 0.020 -0.141 0.179
## 17 lfw 80.421 0.95 -36.714 197.557 1.357 144.054 0.177 0.112 -0.051 0.267
## 18 lfw 16.232 0.95 -110.139 142.603 0.254 145.955 0.800 0.021 -0.140 0.180
## 19 sla -8.929 0.95 -13.382 -4.475 -3.963 143.815 0.000 -0.314 -0.446 -0.160
## 20 sla -9.597 0.95 -14.407 -4.786 -3.943 145.777 0.000 -0.310 -0.442 -0.157
## 21 ldmc 22.906 0.95 17.579 28.233 8.499 145.065 0.000 0.577 0.463 0.663
## 22 ldmc 11.957 0.95 6.268 17.646 4.155 142.344 0.000 0.329 0.176 0.459
## 23 can 179.062 0.95 122.597 235.527 6.268 143.591 0.000 0.464 0.329 0.572
## 24 can -62.650 0.95 -123.768 -1.533 -2.026 145.171 0.045 -0.166 -0.315 -0.004print(st.alltr.cv)## trait estimate stderr bcil bciu df tval pval choenr rcil rciu
## 1 lnc 0.001 0.95 0.000 0.001 2.181 146.886 0.031 0.177 0.017 0.324
## 2 lnc NA NA NA NA NA NA NA NA NA NA
## 3 lcc 0.001 0.95 -0.001 0.003 1.343 143.465 0.181 0.111 -0.052 0.266
## 4 lcc -0.005 0.95 -0.007 -0.003 -5.532 144.973 0.000 -0.418 -0.533 -0.276
## 5 cn 0.032 0.95 0.015 0.049 3.642 144.805 0.000 0.290 0.134 0.424
## 6 cn -0.048 0.95 -0.066 -0.029 -5.087 143.215 0.000 -0.391 -0.512 -0.246
## 7 la 0.180 0.95 0.106 0.255 4.771 143.837 0.000 0.370 0.222 0.493
## 8 la -0.216 0.95 -0.297 -0.135 -5.293 145.831 0.000 -0.401 -0.519 -0.259
## 9 ldw -0.026 0.95 -0.126 0.073 -0.523 140.193 0.602 -0.044 -0.205 0.121
## 10 ldw 0.080 0.95 -0.039 0.198 1.326 137.155 0.187 0.112 -0.055 0.270
## 11 lfw 0.001 0.95 -0.002 0.004 0.649 146.746 0.517 0.053 -0.108 0.211
## 12 lfw NA NA NA NA NA NA NA NA NA NA
## 13 sla -0.291 0.95 -0.380 -0.202 -6.465 143.966 0.000 -0.474 -0.580 -0.341
## 14 sla 0.111 0.95 0.015 0.207 2.281 145.962 0.024 0.186 0.025 0.332
## 15 ldmc -0.249 0.95 -0.321 -0.177 -6.852 143.327 0.000 -0.497 -0.599 -0.367
## 16 ldmc 0.007 0.95 -0.071 0.085 0.180 144.446 0.857 0.015 -0.147 0.175
## 17 can 0.011 0.95 -0.078 0.100 0.245 144.236 0.807 0.020 -0.141 0.181
## 18 can 0.040 0.95 -0.055 0.135 0.832 145.809 0.407 0.069 -0.093 0.225###
interpret_r(st.alltr.cv$choenr, rules = "gignac2016")## [1] "small" NA "small" "large" "moderate" "large" "large"
## [8] "large" "very small" "small" "very small" NA "large" "small"
## [15] "large" "very small" "very small" "very small"
## (Rules: gignac2016)Figures
Make figures
## average ##
for(i in 1:12){
print(plot(Effect(c("years"), mod.alltr.mn[[i]])))
}
dev.off()## null device
## 1####
a1 = ggplot(alltr_mn, aes(x = years,
y = latspr,
colour = site)) +
geom_jitter(width = 0.2) +
geom_smooth(method = lm, formula = y ~ poly(x,2)) +
theme_classic()
a2 = ggplot(alltr_mn, aes(x = years,
y = flost,
colour = site)) +
geom_jitter(width = 0.2) +
geom_smooth(method = lm, formula = y ~ poly(x,2)) +
theme_classic()
a3 = ggplot(alltr_mn, aes(x = years,
y = floper,
colour = site)) +
geom_jitter(width = 0.2) +
geom_smooth(method = lm, formula = y ~ poly(x,2)) +
theme_classic()
a4 = ggplot(alltr_mn, aes(x = years,
y = lnc,
colour = site)) +
geom_jitter(width = 0.2) +
geom_smooth(method = lm, formula = y ~ poly(x,2)) +
theme_classic()
a5 = ggplot(alltr_mn, aes(x = years,
y = lcc,
colour = site)) +
geom_jitter(width = 0.2) +
geom_smooth(method = lm, formula = y ~ poly(x,2)) +
theme_classic()
a6 = ggplot(alltr_mn, aes(x = years,
y = cn,
colour = site)) +
geom_jitter(width = 0.2) +
geom_smooth(method = lm, formula = y ~ poly(x,2)) +
theme_classic()
a7 = ggplot(alltr_mn, aes(x = years,
y = la,
colour = site)) +
geom_jitter(width = 0.2) +
geom_smooth(method = lm, formula = y ~ poly(x,2)) +
theme_classic()
a8 = ggplot(alltr_mn, aes(x = years,
y = ldw,
colour = site)) +
geom_jitter(width = 0.2) +
geom_smooth(method = lm, formula = y ~ poly(x,2)) +
theme_classic()
a9 = ggplot(alltr_mn, aes(x = years,
y = lfw,
colour = site)) +
geom_jitter(width = 0.2) +
geom_smooth(method = lm, formula = y ~ poly(x,2)) +
theme_classic()
a10 = ggplot(alltr_mn, aes(x = years,
y = sla,
colour = site)) +
geom_jitter(width = 0.2) +
geom_smooth(method = lm, formula = y ~ poly(x,2)) +
theme_classic()
a11 = ggplot(alltr_mn, aes(x = years,
y = ldmc,
colour = site)) +
geom_jitter(width = 0.2) +
geom_smooth(method = lm, formula = y ~ poly(x,2)) +
theme_classic()
a12 = ggplot(alltr_mn, aes(x = years,
y = can,
colour = site)) +
geom_jitter(width = 0.2) +
geom_smooth(method = lm, formula = y ~ poly(x,2)) +
theme_classic()
ggarrange(a1, a2, a3, a4, a5, a6, a7, a8, a9, a10, a11, a12,
ncol = 3, nrow = 4)
### cv
for(i in 1:9){
print(plot(Effect(c("years"), mod.alltr.cv[[i]])))
}
dev.off()## null device
## 1a1 = ggplot(alltr_cv, aes(x = years,
y = lnc,
colour = site)) +
geom_jitter(width = 0.2) +
geom_smooth(method = lm, formula = y ~ poly(x,2)) +
theme_classic()
a2 = ggplot(alltr_cv, aes(x = years,
y = lcc,
colour = site)) +
geom_jitter(width = 0.2) +
geom_smooth(method = lm, formula = y ~ poly(x,2)) +
theme_classic()
a3 = ggplot(alltr_cv, aes(x = years,
y = cn,
colour = site)) +
geom_jitter(width = 0.2) +
geom_smooth(method = lm, formula = y ~ poly(x,2)) +
theme_classic()
a4 = ggplot(alltr_cv, aes(x = years,
y = la,
colour = site)) +
geom_jitter(width = 0.2) +
geom_smooth(method = lm, formula = y ~ poly(x,2)) +
theme_classic()
a5 = ggplot(alltr_cv, aes(x = years,
y = ldw,
colour = site)) +
geom_jitter(width = 0.2) +
geom_smooth(method = lm, formula = y ~ poly(x,2)) +
theme_classic()
a6 = ggplot(alltr_cv, aes(x = years,
y = lfw,
colour = site)) +
geom_jitter(width = 0.2) +
geom_smooth(method = lm, formula = y ~ poly(x,2)) +
theme_classic()
a7 = ggplot(alltr_cv, aes(x = years,
y = sla,
colour = site)) +
geom_jitter(width = 0.2) +
geom_smooth(method = lm, formula = y ~ poly(x,2)) +
theme_classic()
a8 = ggplot(alltr_cv, aes(x = years,
y = ldmc,
colour = site)) +
geom_jitter(width = 0.2) +
geom_smooth(method = lm, formula = y ~ poly(x,2)) +
theme_classic()
a9 = ggplot(alltr_cv, aes(x = years,
y = can,
colour = site)) +
geom_jitter(width = 0.2) +
geom_smooth(method = lm, formula = y ~ poly(x,2)) +
theme_classic()
ggarrange(a1, a2, a3, a4, a5, a6, a7, a8, a9,
ncol = 3, nrow = 3)## Warning: Removed 5 rows containing non-finite values (stat_smooth).## Warning: Removed 5 rows containing missing values (geom_point).